Numerical Exercise 3¶

(Deadline: 27/10)¶
Prof. M. Celeste Artale¶

This numerical exercise is based on Valeria Cristiani’s lecture. The goal of this numerical exercise is to familiarize ourselves with the different ways of identifying the components of a galaxy using galaxychop1:

  • JHistogram: Implementation of the dynamic decomposition model of galaxies described by Abadi et al.(2003).

  • JThreshold: Implementation of the dynamic decomposition model of galaxies used by Tissera et al.(2012), Vogelsberger et al.(2014), Marinacci et al.(2014), Park et al.(2019), etc.

  • KMeans: Implementation of Scikit-Learn K-means as a model for the dynamical decomposing of galaxies.

  • GaussianMixture: Implementation of the dynamic decomposition model of galaxies described by Obreja et al.(2018).

  • AutoGaussianMixture: Implementation of the dynamic decomposition model of galaxies described by Du et al.(2019).

Activity¶

1. Install galaxychop. Download the files and scripts available on the NumEx3 folder. The files gal394242.h5, and gal443801.h5 correspond to SubhaloID: 394242 and 443801, respectively, from IllustrisTNG100-1 at z = 0 (snapshot 99)¶

In [1]:
import galaxychop as gchop # Installed and working 
import numpy as np
import matplotlib.pylab as plt
import matplotlib.gridspec as gridspec
import matplotlib as mpl
import pandas as pd


%matplotlib inline

plt.style.use('dark_background')
plt.style.use('seaborn-bright')

COLOR = 'lightgray'

mpl.rcParams['text.color']      = COLOR  # labels, titles, etc colors
mpl.rcParams['axes.labelcolor'] = COLOR  # Numbers in axis colors
mpl.rcParams['xtick.color']     = COLOR  # ticks colors
mpl.rcParams['ytick.color']     = COLOR  # ticks colors
mpl.rc(      'axes',edgecolor   = COLOR) # axis Colors

mpl.rc('font', size=15) # Font Size

mpl.rcParams['figure.figsize'] = [10, 10] # Size of figures
In [2]:
gal  = gchop.read_hdf5("gal394242.h5")
gal2 = gchop.read_hdf5("gal443801.h5")

print(f'SubhaloID 394242: {gal}\nSubhaloID 443801: {gal2}\n')
SubhaloID 394242: <Galaxy stars=37393, dark_matter=155101, gas=80153, potential=True>
SubhaloID 443801: <Galaxy stars=47789, dark_matter=282694, gas=4459, potential=True>

In [3]:
if not gchop.preproc.is_centered_and_aligned(gal):    
    print("Centered & Aligned:", gchop.preproc.is_centered_and_aligned(gal))
    #gal = gchop.preproc.center_and_align(gal)
    gal = gchop.preproc.center_and_align(gal, r_cut=150)
    print("Centered & Aligned:", gchop.preproc.is_centered_and_aligned(gal))
Centered & Aligned: False
Centered & Aligned: True
In [4]:
if not gchop.preproc.is_centered_and_aligned(gal2):    
    print("Centered & Aligned:", gchop.preproc.is_centered_and_aligned(gal2))
    #gal2 = gchop.preproc.center_and_align(gal2)
    gal2 = gchop.preproc.center_and_align(gal2, r_cut=260)
    print("Centered & Aligned:", gchop.preproc.is_centered_and_aligned(gal2))
Centered & Aligned: False
Centered & Aligned: True

2. Using jupyterlab in the IllustrisTNG website, investigate the galaxy properties of the two galaxies: stellar masses, star formation rate, metallicity, gas mass and virial mass of the dark matter halo. Are these galaxies central or satellite galaxies?¶

In the jupyter lab I have run:

########################################

# Packages 

import numpy as np
import matplotlib.pylab as plt
import illustris_python as il
import matplotlib.gridspec as gridspec

########################################

# Getting the data

basePath = '/home/tnguser/sims.TNG/TNG100-1/output/'
snapNum  = 99  # z=0

fields   = ['SubhaloFlag', 'SubhaloHalfmassRadType',
            'SubhaloMassType', 'SubhaloPos', 'SubhaloVel',
            'SubhaloSFRinRad', 'SubhaloSFR','SubhaloStarMetallicity','SubhaloGrNr']

Subhalo  = il.groupcat.loadSubhalos(basePath=basePath,snapNum=snapNum,fields=fields)

FIELDS=['GroupFirstSub','GroupMassType','Group_M_Crit200','Group_R_Crit200']
HHH=il.groupcat.loadHalos(basePath,99,fields=FIELDS)


GFSub = HHH['GroupFirstSub']
Mtype = HHH['GroupMassType']
M200  = HHH['Group_M_Crit200']
R200  = HHH['Group_R_Crit200']


########################################

# Saving the data

gdms=np.array([0,1,4]) # to keep just columns of gas,dm and stars

Data={}
for i in [394242, 443801]:
    a=Subhalo['SubhaloGrNr'][i]
    Data[str(i)]= dict( id          = i, # The subhalo "ID" is the same as the "index"
                        cpos        = Subhalo['SubhaloPos'][i],
                        cvel        = Subhalo['SubhaloVel'][i],
                        rh_type     = Subhalo['SubhaloHalfmassRadType'][i][gdms],
                        M_type      = Subhalo['SubhaloMassType'][i][gdms],
                        sfr_2rh     = Subhalo['SubhaloSFRinRad'][i],     # in 2 * R_half
                        sfr_tot     = Subhalo['SubhaloSFR'][i],          # in 2 * R_half
                        mh          = Subhalo['SubhaloStarMetallicity'][i],
                        Central     = (i in GFSub),
                        M200_parent = M200[a],
                        R200_parent = R200[a])

np.save('Gal_Data.npy', Data)

Now from this file we can print the info

In [5]:
Data=np.load('Gal_Data.npy',allow_pickle=True).item()
In [6]:
pd.DataFrame(Data)
Out[6]:
394242 443801
id 394242 443801
cpos [63794.09, 12026.507, 5045.7017] [15467.963, 63014.26, 73353.84]
cvel [194.05467, 263.4645, 114.92622] [157.85007, -136.74078, -32.544792]
rh_type [60.96873, 65.7992, 5.139353] [109.72671, 74.55964, 1.619801]
M_type [8.251568, 78.41508, 2.5353656] [0.48876184, 142.92282, 3.0524666]
sfr_2rh 0.355641 0.0
sfr_tot 0.93635 0.0
mh 0.01982 0.024343
Central False True
M200_parent 190.210754 139.27948
R200_parent 201.506332 181.628571

3. Using the script Practico.ipynb, visualize the distribution of gas, dm and stellar particles in real space. Discuss the results found for both galaxies¶

In [7]:
x, y, z = gal.stars.x.value , gal.stars.y.value , gal.stars.z.value

r= (x**2+y**2+z**2)**0.5

# print(sum(r>70))

gal.stars.to_dataframe()
Out[7]:
ptype ptypev m x y z vx vy vz softening potential kinetic_energy total_energy Jx Jy Jz
0 stars 0 5.224283e+05 -0.046012 -0.062069 -0.262818 6.362919 22.511274 -5.774254 0.0 -195699.620206 290.293111 -195409.327095 6.274766 -1.937972 -0.640845
1 stars 0 9.745897e+05 -0.155791 -0.125822 -0.183540 21.987698 8.126047 0.272632 0.0 -196176.962277 274.782915 -195902.179362 1.457155 -3.993156 1.500567
2 stars 0 6.935776e+05 -0.250341 -0.162899 -0.206341 -9.982170 7.308657 15.308136 0.0 -195152.120168 193.699612 -194958.420557 -0.985599 5.891984 -3.455738
3 stars 0 1.070959e+06 -0.223984 -0.253224 -0.378401 -9.688729 -5.382013 -1.178408 0.0 -194695.767625 62.113089 -194633.654536 -1.738158 3.402281 -1.247935
4 stars 0 6.013803e+05 -0.068238 -0.198522 -0.285535 17.203435 7.988667 -16.385347 0.0 -195703.942688 314.128285 -195389.814402 5.533893 -6.030286 2.870127
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
37388 stars 0 9.867470e+05 134.949071 -161.114980 -11.970927 -59.939303 57.887495 120.366415 0.0 -26553.893402 10715.877926 -15838.015476 -18699.865471 -15525.806786 -1845.255888
37389 stars 0 8.290853e+05 15.012682 15.493519 -4.697090 -13.645752 262.439571 314.563208 0.0 -101280.155757 84005.373302 -17274.782455 6106.393180 -4658.342086 4151.342534
37390 stars 0 1.004103e+06 17.379644 56.147464 -2.132304 -41.672564 304.593269 11.333777 0.0 -64131.846872 47321.058282 -16810.788589 1285.848336 -108.118416 7633.531240
37391 stars 0 1.526682e+06 269.301518 -82.400695 23.001600 -88.927482 -50.437606 -26.101538 0.0 -19698.409864 5566.669709 -14131.740154 3310.930548 4983.709501 -20910.610228
37392 stars 0 1.356703e+06 13.433154 -8.906366 3.706972 172.602906 380.189578 -150.006406 0.0 -111482.516414 98418.900084 -13063.616330 -73.340046 2654.893305 6644.409980

37393 rows × 16 columns

In [8]:
PLOT=True

LIM=6

if PLOT: # High Resolution
    x, y, z = gal.stars.x.value , gal.stars.y.value , gal.stars.z.value
    
    mask    = (((x**2+y**2+z**2)**0.5)<LIM)
    x, y, z = x[mask], y[mask] , z[mask] 
    
    H, edges=np.histogramdd((x,y,z),bins=50) #bins=(50,45,55)

    x= (edges[0][1:] + edges[0][:-1])/2
    y= (edges[1][1:] + edges[1][:-1])/2
    z= (edges[2][1:] + edges[2][:-1])/2

    x,y,z = np.meshgrid(x,y,z)

    x,y,z=x.ravel(),y.ravel(),z.ravel()

    #------------------------------------------
    x2, y2, z2 = gal2.stars.x.value , gal2.stars.y.value , gal2.stars.z.value
    
    mask       = (((x2**2+y2**2+z2**2)**0.5)<LIM)
    x2, y2, z2 = x2[mask], y2[mask] , z2[mask] 
    

    H2, edges=np.histogramdd((x2,y2,z2),bins=50) #bins=(50,45,55)

    x2= (edges[0][1:] + edges[0][:-1])/2
    y2= (edges[1][1:] + edges[1][:-1])/2
    z2= (edges[2][1:] + edges[2][:-1])/2

    x2,y2,z2 = np.meshgrid(x2,y2,z2)

    x2,y2,z2=x2.ravel(),y2.ravel(),z2.ravel()
In [9]:
if PLOT: # Scatter density version (plotly)
    import plotly.graph_objects as go
    from plotly.subplots import make_subplots

    #------------------
    # Figure

    fig = make_subplots(rows=1, cols=2,
                    specs=[[{'type': 'volume'},{'type': 'volume'}]],
                    print_grid=False,
                    subplot_titles=( ['SubHalo: 394242','SubHalo: 443801']),
                    horizontal_spacing = 0.075
                    )

    #------------------
    # Scatter plot High resolution

    trace1 = go.Volume(name         = '394242',
                       x            = y,
                       y            = x,
                       z            = z,
                       value        = H.ravel(),
                       isomin       = 1,
                       isomax       = H.max(),
                       opacity      = 0.15, # needs to be small to see through all surfaces
                       opacityscale = 'uniform',#""max",
                       surface_count= 25, # needs to be a large number for good volume rendering
                       colorscale   = 'matter',
                       showscale    = False,
#                        colorbar     = dict(thickness    = 20,
#                                            outlinewidth = 0.7,
#                                            title        = "Count particles per bin",
#                                            titleside    = "right",
#                                            tickmode     = "array",
#                                            tickvals     = [], #[0.001,0.999]
#                                            #labelalias   = {1: str(H.max()), 0:str(H.min()) },
#                                            ticks        = "outside",
#                                            x            = -0.2,)  # position of colorbar
                      )
    
    
    #------------------
    # Scatter plot High resolution

    trace2 = go.Volume(name         = '443801',
                       x            = y2,
                       y            = x2,
                       z            = z2,
                       value        = H2.ravel(),
                       isomin       = 1,
                       isomax       = H2.max(),
                       opacity      = 0.15, # needs to be small to see through all surfaces
                       opacityscale = 'uniform',#"max",
                       surface_count= 25, # needs to be a large number for good volume rendering
                       colorscale   = 'matter',
                       colorbar     = dict(thickness    = 20,
                                           outlinewidth = 0.7,
                                           title        = "Count Stars",
                                           titleside    = "right",
                                           tickmode     = "array",
                                           tickvals     = [], #[0.001,0.999]
                                           #labelalias   = {1: str(H.max()), 0:str(H.min()) },
                                           ticks        = "outside",
                                           x            = 1.2,)  # position of colorbar
                      )
    
    fig.add_trace(trace1, # trace= : figure
                  row=1,  # row=   : Position y
                  col=1)  # col=   : Position x
    
    fig.add_trace(trace2, # trace= : figure
                  row=1,  # row=   : Position y
                  col=2)  # col=   : Position x
    
    #------------------
    # Background color (page and axes)

    # https://www.w3schools.com/colors/colors_rgb.asp

    fig['layout']['paper_bgcolor']="rgb(17,17,17)"
        
    set_bgcolor= dict(showbackground  = True,
                      backgroundcolor = "rgb(101, 101, 101)", #"rgb(212, 219, 249)"
                      gridcolor       = "rgb(10, 10, 10)",
                      zeroline        = False)

    fig.update_scenes(xaxis=set_bgcolor,
                      yaxis=set_bgcolor,
                      zaxis=set_bgcolor)
    
    #------------------
    # Font color and size

    fig.update_layout(font_color = 'white', font_size = 16)

    
    fig.update_layout(height=550, width=1000,title_text="Stars in Halos")
    fig.show()
    print()

In [10]:
g = gal.plot(corner=True)

L=250

g.axes[1, 0].set_xlim(-L, L)
g.axes[2, 0].set_xlim(-L, L)
g.axes[2, 1].set_xlim(-L, L)

g.axes[1, 0].set_ylim(-L, L)
g.axes[2, 0].set_ylim(-L, L)
g.axes[2, 1].set_ylim(-L, L)
g.fig.suptitle('SubHalo: 394242')
plt.show()
print()

In [11]:
g = gal2.plot(corner=True)

L=350

g.axes[1, 0].set_xlim(-L, L)
g.axes[2, 0].set_xlim(-L, L)
g.axes[2, 1].set_xlim(-L, L)

g.axes[1, 0].set_ylim(-L, L)
g.axes[2, 0].set_ylim(-L, L)
g.axes[2, 1].set_ylim(-L, L)
g.fig.suptitle('SubHalo: 443801')
plt.show()
print()

The main difference that I can appreciate from these figures between these galaxies are the size, the amount of gas and the shape of the distribution of the stars

4. Following Practico.ipynb, use the JHistogram method for both galaxies and compare the results obtained. Discuss the stellar mass budget in the disk and bulge.¶

In [12]:
decomp  = gchop.models.JThreshold(reassign=True)
In [13]:
components  = decomp.decompose(gal)
In [14]:
print('SubHalo: 394242')
components.describe()
SubHalo: 394242
Out[14]:
Particles Deterministic mass
Size Fraction Size Fraction
Spheroid 12639 0.338013 1.210456e+10 0.323421
Disk 24753 0.661987 2.532209e+10 0.676579
In [15]:
gal.plot.sdyn_pairplot(labels=components, sdyn_kws={'reassign': True}, corner=True)
Out[15]:
<seaborn.axisgrid.PairGrid at 0x141983640>
In [16]:
g = gal.plot.pairplot(labels=components, palette="Pastel2", corner=True)

g.axes[1, 0].set_xlim(-50, 50)
g.axes[2, 0].set_xlim(-50, 50)
g.axes[2, 1].set_xlim(-50, 50)

g.axes[1, 0].set_ylim(-50, 50)
g.axes[2, 0].set_ylim(-50, 50)
g.axes[2, 1].set_ylim(-50, 50)
g.axes[2, 2].set_ylim(-50, 50)
g.fig.suptitle('SubHalo: 394242')
plt.show()
print()

For this galaxy we have almost the same amount of disk than bulge

In [17]:
decomp2     = gchop.models.JThreshold(reassign=True)
components2 = decomp2.decompose(gal2)
print('\nSubHalo: 443801')
components2.describe()
SubHalo: 443801
Out[17]:
Particles Deterministic mass
Size Fraction Size Fraction
Spheroid 39829 0.833434 3.740648e+10 0.830121
Disk 7960 0.166566 7.655026e+09 0.169879
In [18]:
gal2.plot.sdyn_pairplot(labels=components2, sdyn_kws={'reassign': True}, corner=True)
print()

In [19]:
g = gal2.plot.pairplot(labels=components2, palette="Pastel2", corner=True) 

g.axes[1, 0].set_xlim(-50, 50)
g.axes[2, 0].set_xlim(-50, 50)
g.axes[2, 1].set_xlim(-50, 50)

g.axes[1, 0].set_ylim(-50, 50)
g.axes[2, 0].set_ylim(-50, 50)
g.axes[2, 1].set_ylim(-50, 50)
g.axes[2, 2].set_ylim(-50, 50)
g.fig.suptitle('SubHalo: 443801')
plt.show()
print()

while for this galaxy we have complete dominated by spheroidal component

5. Investigate the distribution of the stellar ages of each component, and make a plot of the circularity parameter vs. stellar age for the particles. You will need to use the additional files age_gal443801.pkl and age_gal394242.pkl For reading the files you can use:¶

import pandas as pd
df = pd.read_pickle("age_gal394242.pkl")
In [20]:
import pandas as pd
df = pd.read_pickle("age_gal394242.pkl")
age1=df.to_numpy().T[0]

df = pd.read_pickle("age_gal443801.pkl")
age2=df.to_numpy().T[0]

del df
In [21]:
eps1=gal.stellar_dynamics().eps
eps2=gal2.stellar_dynamics().eps
In [22]:
from sklearn.neighbors import KernelDensity # I used a KDE to improve the view of the distirbution

A_plot = np.linspace(-0.2, 14 , 2000)[:, np.newaxis]
e_plot = np.linspace(-1.1, 1.1, 2000)[:, np.newaxis]


A1      = age1[~np.isnan(eps1)][:, np.newaxis]
kde     = KernelDensity(kernel="cosine", bandwidth=0.033).fit(A1)
A1_dens = np.exp(kde.score_samples(A_plot))

A2      = age2[~np.isnan(eps2)][:, np.newaxis]
kde     = KernelDensity(kernel="cosine", bandwidth=0.033).fit(A2)
A2_dens = np.exp(kde.score_samples(A_plot))

e1      = eps1[~np.isnan(eps1)][:, np.newaxis]
kde     = KernelDensity(kernel="cosine", bandwidth=0.033).fit(e1)
e1_dens = np.exp(kde.score_samples(e_plot))

e2      = eps2[~np.isnan(eps2)][:, np.newaxis]
kde     = KernelDensity(kernel="cosine", bandwidth=0.033).fit(e2)
e2_dens = np.exp(kde.score_samples(e_plot))
In [23]:
fig= plt.figure(figsize=(20,10))
gs = gridspec.GridSpec( 2,5,                             # Number of axis y,x
                        height_ratios = [0.3,1],             # relatives ratios of heigh
                        width_ratios  = [0.3,1,0.15,1,0.3], # relatives ratio of with
                        left  = 0.1,      # Space to the edge of left   from the nearest axis
                        right = 0.95,     # Space to the edge of right  from the nearest axis
                        bottom= 0.15,     # Space to the edge of bottom from the nearest axis
                        top   = 0.97,     # Space to the edge of top    from the nearest axis
                        wspace= 0.02,      # Space horizontal between each of the axis
                        hspace= 0.02)      # Space vertical   between each of the axis

C1   = fig.add_subplot(gs[0,0])
C2   = fig.add_subplot(gs[0,4])

ax1  = fig.add_subplot(gs[1,1])
ax3  = fig.add_subplot(gs[0,1],sharex=ax1)
ax5  = fig.add_subplot(gs[1,0],sharey=ax1)

ax2  = fig.add_subplot(gs[1,3],sharey=ax1,sharex=ax1)
ax4  = fig.add_subplot(gs[0,3],sharex=ax1)
ax6  = fig.add_subplot(gs[1,4],sharey=ax1)



#------------------------------------------
# Corners saying density
C1.text(0.5,0.5,'Density',verticalalignment='center',horizontalalignment='center',rotation=45)
plt.setp(C1.spines.values(), visible=False)
C1.minorticks_off()
C1.set_xticks([])
C1.set_yticks([])


C2.text(0.5,0.5,'Density',verticalalignment='center',horizontalalignment='center',rotation=-45)
plt.setp(C2.spines.values(), visible=False)
C2.minorticks_off()
C2.set_xticks([])
C2.set_yticks([])

#------------------------------------------
# Gal 394242

ax1.scatter(age1,eps1,marker='.',alpha=0.1)

ax3.plot(A_plot,A1_dens)
ax3.set_title('Gal 394242')

ax5.plot(e1_dens,e_plot)
ax5.invert_xaxis()

#------------------------------------------
# Gal 433801

ax2.scatter(age2,eps2,marker='.',alpha=0.07)

ax4.plot(A_plot,A2_dens)
ax4.set_title('Gal 433801')

ax6.plot(e2_dens,e_plot)

#------------------------------------------
# Erasing extra label in axis

plt.setp(ax1.get_yticklabels(), visible=False)
plt.setp(ax3.get_xticklabels(), visible=False)

plt.setp(ax4.get_xticklabels(), visible=False)
plt.setp(ax6.get_yticklabels(), visible=False)

#------------------------------------------
# Axes labels
fig.supylabel('$J_z$ / $J_{circ}$',fontsize=22)
fig.supxlabel('Age (Gyr)',fontsize=22)
plt.show()
print()

6. Extend your analysis of the stellar components by using one of the other alternative methods available in galaxychop. Examples:¶


#Initialize the decomposer KMeans:
decomposer = gchop.models.KMeans(n_components=2, reassign=True) #OR Initialize the decomposer GaussianMixture:
decomposer = gchop.models.GaussianMixture(n_components=2, reassign=True)
components = decomposer.decompose(gal).to_dataframe()
In [24]:
decomposer1 = gchop.models.GaussianMixture(n_components=2, reassign=True)
components1 = decomposer1.decompose(gal)

decomposer2 = gchop.models.GaussianMixture(n_components=2, reassign=True)
components22= decomposer2.decompose(gal2)
In [25]:
gal.plot.sdyn_pairplot(labels=components1, sdyn_kws={'reassign': True}, corner=True)
print()

In [26]:
g = gal.plot.pairplot(labels=components1, palette="Pastel2", corner=True)

g.axes[1, 0].set_xlim(-50, 50)
g.axes[2, 0].set_xlim(-50, 50)
g.axes[2, 1].set_xlim(-50, 50)

g.axes[1, 0].set_ylim(-50, 50)
g.axes[2, 0].set_ylim(-50, 50)
g.axes[2, 1].set_ylim(-50, 50)
g.axes[2, 2].set_ylim(-50, 50)
g.fig.suptitle('SubHalo: 394242')
plt.show()
print()

In [27]:
gal2.plot.sdyn_pairplot(labels=components22, sdyn_kws={'reassign': True}, corner=True)
print()

In [28]:
g = gal2.plot.pairplot(labels=components22, palette="Pastel2", corner=True) 

g.axes[1, 0].set_xlim(-50, 50)
g.axes[2, 0].set_xlim(-50, 50)
g.axes[2, 1].set_xlim(-50, 50)

g.axes[1, 0].set_ylim(-50, 50)
g.axes[2, 0].set_ylim(-50, 50)
g.axes[2, 1].set_ylim(-50, 50)
g.axes[2, 2].set_ylim(-50, 50)
g.fig.suptitle('SubHalo: 443801')
plt.show()
print()

In [29]:
print("\n\n\n\nBy Ian Baeza")



By Ian Baeza